rm(list=ls())
pacman::p_load(foreign, ggplot2, ggpubr, plm, reshape2, countrycode, sandwich, lmtest, MASS, 
               rworldmap, RColorBrewer, states, mice, VIM, stargazer, margins, clusterSEs, lme4, optimx,
               coefplot, gridExtra, stringr, xtable, plyr, dplyr, viridis)

load("covid_wave1_merged.RData")

##### DESCRIPTIVE STATISTICS AND FIGURES #####

# Table 1
table(df$covid_contract_dum,df$covid_job_dum)
prop.table(table(df$covid_contract_dum,df$covid_job_dum))

# Table 2 
table(df$covid_contract_dum,df$partyid2)
prop.table(table(df$covid_contract_dum,df$partyid2))
table(df$covid_job_dum,df$partyid2)
prop.table(table(df$covid_job_dum,df$partyid2))

# Figure 2
ggplot(df) + geom_bar(aes(x=covid_help)) + geom_text(aes(x=covid_help,label=..count..),stat="count",vjust=-1)+
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank(),axis.text.x = element_blank()) +
  ggtitle("On a scale of 1 to 10, should the top priority of the U.S. be to...")+labs(fill="") + annotate("text",x=.5,y=-50,hjust=0,label = "1 - Respond to COVID-19\n at home") + ylim(c(-50,530))+
  annotate("text",x=10.5,y=-50,hjust=1,label = "10 - Help developing countries\n respond to COVID-19")
ggsave("covid_help.pdf", width =9)

# Figure 4
a <- ggplot(subset(df,!is.na(partyid2))) + geom_bar(aes(x=partyid2,fill=covid_help), position = "fill") +
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank()) +
  labs(fill="")
b <- ggplot(subset(df,!is.na(covid_contract))) + geom_bar(aes(x=covid_contract,fill=covid_help), position = "fill") +
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank()) +
  labs(fill="")
p <- ggarrange(a,b,nrow=1,common.legend = TRUE, legend="left")
annotate_figure(p,top="Should the top priority of the U.S. be to...")
ggsave("covid_help_determinants.pdf",width=10)

# Figure 1
ggplot(df) + geom_bar(aes(x=covid_contrib)) + geom_text(aes(x=covid_contrib,label=..count..),stat="count",vjust=-1)+
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank()) +
  ggtitle("Should the U.S. contribute money to international organizations like the WHO?")+labs(fill="")+ylim(-30,900)
ggsave("covid_contrib.pdf",width = 10)

# Figure 3 
a <- ggplot(subset(df,!is.na(partyid2))) + geom_bar(aes(x=partyid2,fill=covid_contrib), position = "fill") +
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank()) +
  labs(fill="")
b <- ggplot(subset(df,!is.na(covid_contract))) + geom_bar(aes(x=covid_contract,fill=covid_contrib), position = "fill") +
  theme_bw() + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title.x = element_blank(),axis.title.y=element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.ticks.x=element_blank()) +
  labs(fill="")
p <- ggarrange(a,b,nrow=1,common.legend = TRUE, legend="right")
annotate_figure(p,top="Should the U.S. contribute money to international organizations like the WHO?")
ggsave("covid_contrib_determinants.pdf",width=8)

# Table 5
m1 <- lm(data = df, formula = covid_contract_dum ~ republican)
m2 <- lm(data = df, formula = covid_contract_dum ~ republican+income_num+female+edu_num+age+nonwhite+state)
m3 <- lm(data = df, formula = covid_contract_dum ~ nyt_log_cases_apr_pc+th_per_gop)
se3 <- coeftest(m3, vcov = vcovCL(m3, cluster = df$fips))
m4 <- lm(data = df, formula = covid_contract_dum ~ republican+income_num+female+edu_num+age+nonwhite+nyt_log_cases_apr_pc+th_per_gop)
se4 <- coeftest(m4, vcov = vcovCL(m4, cluster = df$fips))
m5 <- lm(data = df, formula = covid_contract_dum ~ republican+income_num+female+edu_num+age+nonwhite+nyt_log_cases_apr_pc+th_per_gop*republican)
se5 <- coeftest(m5, vcov = vcovCL(m5, cluster = df$fips))
m6 <- lm(data = df, formula = covid_job_dum ~ republican)
m7 <- lm(data = df, formula = covid_job_dum ~ republican+income_num+female+edu_num+age+nonwhite+state)
m8 <- lm(data = df, formula = covid_job_dum ~ nyt_log_cases_apr_pc+th_per_gop)
se8 <- coeftest(m8, vcov = vcovCL(m8, cluster = df$fips))
m9 <- lm(data = df, formula = covid_job_dum ~ republican+income_num+female+edu_num+age+nonwhite+nyt_log_cases_apr_pc+th_per_gop)
se9 <- coeftest(m9, vcov = vcovCL(m9, cluster = df$fips))
m10 <- lm(data = df, formula = covid_job_dum ~ republican+income_num+female+edu_num+age+nonwhite+nyt_log_cases_apr_pc+th_per_gop*republican)
se10 <- coeftest(m10, vcov = vcovCL(m10, cluster = df$fips))
cat(stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,
       se = list(summary(m1)$coefficients[,2],summary(m2)$coefficients[,2],se3[,2],se4[,2],se5[,2],summary(m6)$coefficients[,2],summary(m7)$coefficients[,2],se8[,2],se9[,2],se10[,2]),
       dep.var.labels = c("Contracted (0-1)","Lost Job (0-1)"),
       covariate.labels = c("Republican (0-1)",
                            "Income",
                            "Female",
                            "Edu",
                            "Age",
                            "Nonwhite",
                            "County Cases (log pc)",
                            "County Trump VS",
                            "Rep*County Trump VS"),
       title = "Correlation between partisanship and exposure to COVID-19",
       label = "covid_republican",
       font.size = "small",no.space = TRUE,
       add.lines = list(c("Fixed Effects?","none","state","none","none","none","none","state","none","none","none"),
                        c("Clustered S.E.s?","none","none","county","county","county","none","none","county","county","county")),
       omit = "state", omit.stat = c("ll","rsq","ser","f")),file="covid_republican.tex", sep="\n")

##### ANALYSIS #####

m1 <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)
m2 <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)
m3 <- lm(data = df, formula = covid_contrib_num~partyid2+covid_job_dum+income_num+female+edu_num+age+nonwhite+state)
m4 <- lm(data = df, formula = covid_contrib_num~partyid2+covid_job_dum+partyid2*covid_job_dum+income_num+female+edu_num+age+nonwhite+state)
m5 <- lm(data = df, formula = covid_contrib_num~partyid2+nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state)
se5 <- coeftest(m5, vcov = vcovCL(m5, cluster = df$fips))
m6 <- lm(data = df, formula = covid_contrib_num~partyid2*nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state)
se6 <- coeftest(m6, vcov = vcovCL(m6, cluster = df$fips))

m7 <- lm(data = df, formula = covid_help_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)
m8 <- lm(data = df, formula = covid_help_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)
m9 <- lm(data = df, formula = covid_help_num~partyid2+covid_job_dum+income_num+female+edu_num+age+nonwhite+state)
m10 <- lm(data = df, formula = covid_help_num~partyid2+covid_job_dum+partyid2*covid_job_dum+income_num+female+edu_num+age+nonwhite+state)
m11 <- lm(data = df, formula = covid_help_num~partyid2+nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state)
se11 <- coeftest(m11, vcov = vcovCL(m11, cluster = df$fips))
m12 <- lm(data = df, formula = covid_help_num~partyid2*nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state)
se12 <- coeftest(m12, vcov = vcovCL(m12, cluster = df$fips))

# Table 3
cat(stargazer(m1,m2,m3,m4,m5,m6,
              se = list(summary(m1)$coefficients[,2],summary(m2)$coefficients[,2],summary(m3)$coefficients[,2],summary(m4)$coefficients[,2],se5[,2],se6[,2]),
              dep.var.labels = "IO Contributions (1-5)",
              order = c("partyid2","covid_contract_dum","covid_job_dum","nyt_log_cases_apr_pc"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Party:Ind.*Lost Job","Party:Rep.*Lost Job","Party:Ind.*County Cases","Party:Rep.*County Cases",
                                   "Contracted (0-1)",
                                   "Lost Job (0-1)","County Cases (log pc)",
                                   "Income","Female","Edu","Age","Nonwhite"),
              title = "Predictors of attitudes toward IO contributions",
              label = "io_contributions",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","state","state","state","state","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","county","county")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="io_contributions.tex", sep="\n")

# Table 4
cat(stargazer(m7,m8,m9,m10,m11,m12,
              se = list(summary(m7)$coefficients[,2],summary(m8)$coefficients[,2],summary(m9)$coefficients[,2],summary(m10)$coefficients[,2],se11[,2],se12[,2]),
              dep.var.labels = c("Foreign Aid (1-10)"),
              order = c("partyid2","covid_contract_dum","covid_job_dum"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Party:Ind.*Lost Job","Party:Rep.*Lost Job","Party:Ind.*County Cases","Party:Rep.*County Cases",
                                   "Contracted (0-1)",
                                   "Lost Job (0-1)","County Cases (log pc)",
                                   "Income","Female","Edu","Age","Nonwhite"),
              title = "Predictors of attitudes toward foreign aid",
              label = "aid_developing",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","state","state","state","state","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","county","county")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="aid_developing.tex", sep="\n")

# Robustness check: introduce controls gradually

m1r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum)
m2r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum)
m3r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite)
m4r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite)
m5r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)
m6r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state)

# Table 6 (Appendix)
cat(stargazer(m1r,m2r,m3r,m4r,m5r,m6r,
              dep.var.labels = "IO Contributions (1-5)",
              order = c("partyid2","covid_contract_dum"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Contracted (0-1)",
                                   "Income","Female","Edu","Age","Nonwhite"),
              title = "Predictors of attitudes toward IO contributions (controls introduced gradually)",
              label = "io_contributions_gradualcontrols",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","none","none","none","none","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","none","none")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="io_contributions_gradualcontrols.tex", sep="\n")

# Robustness check: are partisan differences driven by political geography?

m1r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m2r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m3r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m4r <- lm(data = df, formula = covid_contrib_num~partyid2+covid_job_dum+partyid2*covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m5r <- lm(data = df, formula = covid_contrib_num~partyid2+nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
se5r <- coeftest(m5r, vcov = vcovCL(m5r, cluster = df$fips))
m6r <- lm(data = df, formula = covid_contrib_num~partyid2*nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
se6r <- coeftest(m6r, vcov = vcovCL(m6r, cluster = df$fips))

m7r <- lm(data = df, formula = covid_help_num~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m8r <- lm(data = df, formula = covid_help_num~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m9r <- lm(data = df, formula = covid_help_num~partyid2+covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m10r <- lm(data = df, formula = covid_help_num~partyid2+covid_job_dum+partyid2*covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m11r <- lm(data = df, formula = covid_help_num~partyid2+nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
se11r <- coeftest(m11r, vcov = vcovCL(m11r, cluster = df$fips))
m12r <- lm(data = df, formula = covid_help_num~partyid2*nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
se12r <- coeftest(m12r, vcov = vcovCL(m12r, cluster = df$fips))

# Table 7 (Appendix)
cat(stargazer(m1r,m2r,m3r,m4r,m5r,m6r,
              se = list(summary(m1r)$coefficients[,2],summary(m2r)$coefficients[,2],summary(m3r)$coefficients[,2],summary(m4r)$coefficients[,2],se5r[,2],se6r[,2]),
              dep.var.labels = "IO Contributions (1-5)",
              order = c("partyid2","covid_contract_dum","covid_job_dum","nyt_log_cases_apr_pc"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Party:Ind.*Lost Job","Party:Rep.*Lost Job","Party:Ind.*County Cases","Party:Rep.*County Cases",
                                   "Party:Ind.*County Trump VS","Party:Rep.*County Trump VS",
                                   "Contracted (0-1)",
                                   "Lost Job (0-1)","County Cases (log pc)",
                                   "Income","Female","Edu","Age","Nonwhite",
                                   "County Trump VS"),
              title = "Robustness: Partisan differences in support for IOs not driven by political geography",
              label = "io_contributions_r",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","state","state","state","state","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","county","county")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="io_contributions_r.tex", sep="\n")

# Table 8 (Appendix)
cat(stargazer(m7r,m8r,m9r,m10r,m11r,m12r,
              se = list(summary(m7r)$coefficients[,2],summary(m8r)$coefficients[,2],summary(m9r)$coefficients[,2],summary(m10r)$coefficients[,2],se11r[,2],se12r[,2]),
              dep.var.labels = c("Foreign Aid (1-10)"),
              order = c("partyid2","covid_contract_dum","covid_job_dum"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Party:Ind.*Lost Job","Party:Rep.*Lost Job","Party:Ind.*County Cases","Party:Rep.*County Cases",
                                   "Party:Ind.*County Trump VS","Party:Rep.*County Trump VS",
                                   "Contracted (0-1)",
                                   "Lost Job (0-1)","County Cases (log pc)",
                                   "Income","Female","Edu","Age","Nonwhite",
                                   "County Trump VS"),
              title = "Robustness: Foreign aid results not driven by political geography",
              label = "aid_developing_r",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","state","state","state","state","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","county","county")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="aid_developing_r.tex", sep="\n")

# Robustness: operationalize DV as a dummy variable to mitigate ceiling effects

m1r <- lm(data = df, formula = covid_contrib_dum~partyid2+covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m2r <- lm(data = df, formula = covid_contrib_dum~partyid2+covid_contract_dum+partyid2*covid_contract_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m3r <- lm(data = df, formula = covid_contrib_dum~partyid2+covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
m4r <- lm(data = df, formula = covid_contrib_dum~partyid2+covid_job_dum+partyid2*covid_job_dum+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
m5r <- lm(data = df, formula = covid_contrib_dum~partyid2+nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop)
se5r <- coeftest(m5r, vcov = vcovCL(m5r, cluster = df$fips))
m6r <- lm(data = df, formula = covid_contrib_dum~partyid2*nyt_log_cases_apr_pc+income_num+female+edu_num+age+nonwhite+state+th_per_gop*partyid2)
se6r <- coeftest(m6r, vcov = vcovCL(m6r, cluster = df$fips))

# Table 9 (Appendix)
cat(stargazer(m1r,m2r,m3r,m4r,m5r,m6r,
              se = list(summary(m1r)$coefficients[,2],summary(m2r)$coefficients[,2],summary(m3r)$coefficients[,2],summary(m4r)$coefficients[,2],se5r[,2],se6r[,2]),
              dep.var.labels = "IO Contributions (0-1)",
              order = c("partyid2","covid_contract_dum","covid_job_dum","nyt_log_cases_apr_pc"),
              covariate.labels = c("Party:Ind.","Party:Rep.","Party:Ind.*Contracted","Party:Rep.*Contracted",
                                   "Party:Ind.*Lost Job","Party:Rep.*Lost Job","Party:Ind.*County Cases","Party:Rep.*County Cases",
                                   "Party:Ind.*County Trump VS","Party:Rep.*County Trump VS",
                                   "Contracted (0-1)",
                                   "Lost Job (0-1)","County Cases (log pc)",
                                   "Income","Female","Edu","Age","Nonwhite",
                                   "County Trump VS"),
              title = "Robustness: Partisan differences in support for IOs not driven by ceiling effects",
              label = "io_contributions_r2",
              font.size = "small",no.space = TRUE,
              add.lines = list(c("Fixed Effects?","state","state","state","state","state","state"),
                               c("Clustered S.E.s?","none","none","none","none","county","county")),
              omit = c("Constant","state"), omit.stat = c("ll","rsq","ser","f")),file="io_contributions_r2.tex", sep="\n")
